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Abstract 

We apply an unsteady Reynolds -averaged Navier-Stokes (URANS) solver for unstructured grids to unsteady flows on 
moving and stationary grids. Example problems considered are relevant to active flow control and stability and control. Com- 
putational results are presented using the Spalart-Allmaras turbulence model and are compared to experimental data. The effect 
of grid and time-step refinement are examined. 


Introduction 

Unsteady flows are being studied with increasing interest both computationally and experimentally. There are many reasons 
for this, but they certainly include the growing realization that unsteady flows exhibit properties that can be exploited for control 
of advanced flight vehicles. In addition, there is simply the fact that all flight vehicles must maneuver at one point or another 
and as a result experience unsteady conditions. NASA has devoted increasing resources to the study of unsteady flows in a 
number of areas. Motivating the current work are NASA research programs focusing on synthetic jets and predictive methods 
for obtaining stability and control (S&C) derivatives. 

In the area of flow control, Thomas, Choudhari, and Joslin 1 have identified several active mechanisms, including synthetic 
jets, as feasible for improving performance and controllability. Potential benefits from active control include reducing drag, 
increasing lift, and enhancing control forces. For non-conventional configurations, such as the blended wing body (BWB), 
active flow control may permit the use of thicker wing sections without adversely affecting performance. 

In the area of stability and control, the emergence of novel aircraft configurations is placing an increased need for predictive 
methods for S&C derivatives beyond traditional methods. Indeed, the classic view of dynamic system linearization about a set 
point with a resulting set of derivatives may be insufficient to deliver the required information for either the design of these new 
configurations or flight simulators for training. 

In an era of increasing costs of experimental testing at a decreasing number of available facilities, there is a growing empha- 
sis on validation of new concepts via computational fluid dynamics (CFD). There have been numerous workshops conducted to 
evaluate the CFD state-of-the-art for predicting flows of interest to the aerodynamics community. A number of workshops (for 
example, the AIAA First 2 and Second 3 CFD Drag Prediction Workshops) have focused on the ability of CFD to predict steady- 
state performance metrics at cruise conditions. More recently, the CFDVAF2004 workshop, 4 had as its primary objective an 
assessment of the the state of the art for measuring and computing aerodynamic flows in the presence of synthetic jets. This 
workshop was comprised of three configurations deemed to be of interest for synthetic jet evaluation. This paper will present 
a CFD simulation of one of the three configurations, designated by the workshop as Case 2. 4 This case represents an isolated 
synthetic jet formed by a single diaphragm-mounted piston actuator exhausting into a nominally uniform cross-flow. Multiple 
measurement techniques, including PIV, FDV, and hot-wire probes were used to generate a large body of experimental data for 
this configuration. References 4 and 5 describe the details of the experimental setup and geometric configuration. 

In September, 2003, NASA Fangley Research Center hosted the first of several planned COMSAC symposia. 6 The inau- 
gural symposium sought to bring together researchers from the stability and control S&C community (typically non experts 
in CFD) and the CFD community (typically non experts in S&C) to discuss issues that currently limit the use of CFD for 
S&C applications and to formulate plans for expanded use of CFD in this area. Subsequent to the symposium, NASA Fangley 
Research Center embarked on a multi-year project to evaluate state-of-the-art CFD codes for calculation of both static and dy- 
namic stability and control characteristics of a Blended Wing Body (BWB) configuration. This paper highlights the application 
of one of the CFD codes under evaluation to problems of forced pitching oscillation, including the BWB. 

Governing Equations 

The unsteady, compressible, Navier-Stokes equations may be written for a moving or stationary control volume in integral 
conservation form as 
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where V is the control volume, bounded by control surface dV, with local control volume face speed |W|. The vector q 
represents the conserved variables; the vectors F* and F v represent the convective and diffusive fluxes of the conserved vari- 
ables, respectively. For the case of a moving control volume, the convective fluxes must account for the relative reduction or 
enhancement of the flux through the control surface owing to the local control surface velocity. Compared with a stationary 
control volume with flux vector F, F* = F — qW. Introducing a volume average of q 


the conservation equations are 
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The components of the stress tensor (with the Stokes hypothesis) and heat flux vector are given by 
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where ji t is a turbulent eddy viscosity determined via a suitable turbulence model. The mean flow equations are closed by the 
equation of state for a perfect gas 
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where p is static pressure and E the total energy per unit volume. The laminar viscosity is obtained via Sutherland’s law 


P = 


1 + C* y3/2 

T + C* 


where C* = 198.6/7’,,,. The governing equations above have been made dimensionless using the free-stream density 
(where ~ denotes a dimensional quantity), speed of sound molecular viscosity , thermal conductivity , and a scaled 
reference length L re f such that the Reynolds number appearing in the stress tensor and heat flux vector is a Reynolds number 
per unit length in the mesh. 


Solution Algorithm 

The unstructured mesh solver used for this study is FUN3D. 7 Within the code, Eq. 1 is discretized over the median dual 
volume surrounding each mesh point, balancing the time rate of change of the averaged conserved variables in each dual volume 
with the flux of mass, momentum, and energy through the instantaneous surface of the control volume. 

Spatial Discretization 

Details of the spatial discretization have been presented in Ref. 7, and will only be summarized here. Within FUN3D, the 
convective fluxes are computed with a flux splitting scheme, and for second order accuracy the values at dual-cell interfaces are 
reconstructed using gradients at mesh nodes computed using a least-squares technique. Limiting of the reconstructed values 
may be employed for flows with strong shocks. For all results presented in this paper, the convective flux scheme used is 
Roe’s flux difference splitting 8 with second order unlimited reconstruction. For tetrahedral meshes, the full viscous fluxes are 
discretized using a finite-volume formulation in which the required velocity gradients on the dual faces are computed using the 
Green-Gauss theorem. On tetrahedral meshes this is equivalent to a Galerkin type approximation. For non-tetrahedral meshes, 
the same Green-Gauss approach can lead to odd-even decoupling. A pure edge-based approach can be used to circumvent 
the odd-even decoupling issue, but yields only approximate viscous terms. Thus for non-tetrahedral meshes, the edge-based 
gradients are combined with Green-Gauss gradients, which improves the h-ellipticity of the operator, and allows the complete 
viscous stresses to be evaluated. For turbulent flows, both the one-equation model of Spalart and Allmaras 9 (SA) and the two- 
equation SST model of Menter 10 are available. The SA model may be solved loosely coupled to the mean-flow equations or 
tightly coupled to the mean-flow equations. For all results presented in this paper, the one equation SA model is employed, 
solved in a loosely coupled fashion. 

Few modifications are needed to the spatial discretization summarized above to handle moving meshes. Specifically, the 
dual face speeds need to be included in the convective flux routines (including boundary routines); the diffusive flux routines 
remain unaltered. For viscous flows, the no-slip condition at the wall is satisfied by setting the fluid velocity equal to the velocity 
of the wall itself; this in turn adds a kinetic energy term to the boundary condition for the energy equation at the wall that does 
not normally appear for stationary meshes. 

Temporal Discretization 

In this paper we restrict ourselves to the class of problems in which the grid either moves as a rigid body or is stationary. 
Thus the control volumes are invariant in time and Eq. 1 may be written as 




Evaluating this equation at time level n+1, and writing the time derivative as a series expansion of successive levels backward 
in time gives 
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The sequence { </„ } governs the accuracy of the temporal discretization; several choices are described below. At this stage the 
right-hand side of Eq. 2 could be linearized about time level n+1, resulting in an implicit scheme for AQ" = Q n+ I — Q". 
However, the linearization introduces an additional error into the result that is time-step dependent. Furthermore, R(Q) is 
difficult to linearize exactly, so that in practice approximate linearizations are often used, introducing another level of error. 

To mitigate these linearization errors, as in Reference 11a pseudo-time term may be added to Eq. 2 
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where the pseudo time is denoted by r and 0 < r < oo during each physical time step. Provided that dQ/dr vanishes for large 
r, during each time step, Eq. 2 is recovered. Discretizing this pseudo-time term with a first order backward difference about 
pseudo-time level m+1, and noting that Q m+1 — > Q” +1 as m — > oo gives 
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where the consistency requirement S <fi n - 0 has been invoked. Equation 3 is the means by which the solution is advanced in 
time in FUN3D; with the physical time step At set to 00 , Eq. 3 reverts to the standard “steady state” scheme described in Ref. 
7. At each subiteration m, the linear system represented by Eq. 3 is iteratively solved using a user-specified number of point 
Jacobi or point Gauss-Seidel sweeps - typically the latter as it converges the linear system approximately twice as fast as point 
Jacobi. More recently a line solver has been added, 12 allowing points near solid walls to be solved in a more closely coupled 
fashion. Within the non-linear iteration process between time steps, the equations are advanced in pseudo time with local time 
stepping (spatially varying time step At based on a constant CFL number) to accelerate the solution to a steady state in pseudo 
time, with the physical time step At being held constant over the entire mesh. The CFL number may be ramped during the 
subiterations in order to accelerate the convergence in pseudo time. To reduce computational time, the evaluation of the flux 
jacobian <9R m /<9Q may be periodically frozen during pseudo-time stepping. 

Returning to the choice of { <j> } „, if one retains terms to 0 n _ 2 , various backward-difference schemes up to third order in 
time may be constructed, as summarized in the following table 
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Table 1 Coefficients for BDF schemes up to third order 


The coefficients for the first-, second- and third-order schemes listed above are well known, with the first- and second-order 
schemes being stable for any time step. The third-order scheme, while usually found to be stable in practice, is not guaranteed 
so. The coefficients for the scheme labeled BDF 2 opt are from Ref. 13, and are a linear combination of second and third order 
coefficients that lead to a scheme that is stable for any time step, with an order of accuracy that is in between second order and 
third order. 
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Temporal Error Control 

The subiterative scheme outlined above has been implemented in a number of flow solvers. One issue that consistently 
arises is the choice of the number of subiterations to use. Converging the subiteration residual at each time step (as measured 
by the complete right-hand side of Eq. 3) to machine zero is generally prohibitive in terms of computational cost. So the issue 
becomes a balance between computational cost and performing just enough subiterations to “get the right answer”. For the BDF 
schemes, we can estimate the temporal error incurred at each time step by examining the difference in residual contribution 
with two different levels of approximations of time derivatives. 14 The temporal error norm can be used as an exit criteria for 
terminating the subiteration loop of the dual time stepping process. Basically, we terminate this subiteration process when the 
residuals (algebraic errors) drop below a specified fraction of the temporal error norm. For all cases presented in this paper, the 
specified fraction was 0.1, and was chosen based on prior experience. Such a strategy results in uniform temporal accuracy for 
all time steps, and eliminates the guess work for selecting iteration count or arbitrary level of residual reduction. 

Mesh Movement 

The initial implementation of time-varying meshes within the FUN3D code were restricted to translation or rotation as a 
solid body. Solid-body translation is trivial, with cell face speeds directly following from the specified translation velocity. 
Rigid rotation is also straightforward, though some details are warranted. 

Consider a rotation, parallel to the y-axis, of a vector A lying in the x-z plane, through an angle A 0, to a new position where 
the vector is denoted by A'. One can show that 


A x > = A x cosA8 + A z sinA9 (4) 

A z i = A z cosA6 — A x sinA6 

Transformations for rotation about the x or z axes are obtained by a cyclic permutation of the coordinate directions in Eq. 4. 
By differentiating Eq. 4 with respect to time, one can show that the corresponding velocity components are 

Ax' = A z ,9 (5) 

A z i = —A x '6 

When applied to a position vector, Eq. 4 may be used to obtain expressions for the new mesh point locations and velocities. For 
example, for rotation parallel to the y-axis, A = r = (x — Xo)i + (z — 2 o)k, with (xo, Zo) the rotation point, gives the required 
expressions for the new position ( x 1 — Xq, z' — zo). 

Because the dual-cell volumes are fixed in time, they may be precomputed and stored. Likewise, dual-face areas do not 
change with time. With rotational motion, dual-face normals do change direction with each time step, but are easily recalculated 
using Eq. 4. The dual-face speeds, needed to evaluate the fluxes via Eq. 1, are calculated from the mesh point velocities obtained 
via Eq. 5 with the same Green-Gauss approach used to evaluate the dual-cell normals for static meshes. 

Boundary Conditions 

The boundary conditions required for solving the Navier-Stokes equations, such as the no-slip, adiabatic wall, far-field 
and extrapolation conditions, are described in Ref. 7, with the above-mentioned extensions of the no-slip condition and the 
convective-flux evaluation for moving grids. 

For the synthetic jet case described below, a special boundary condition was added to simulate the effect of an oscillating 
piston. This is done with a periodic velocity transpiration condition imposed at the mean piston surface. Further details of this 
transpiration condition are given in the section describing the synthetic jet results. 

Pitching Airfoil 

As a first test of the new moving mesh capability within FUN3D, a pitching 0012 airfoil is considered. The flow and 
geometry motion conditions are chosen to represent AGARD Case 3, 15 for which experimental data are available. The airfoil 
is pitched about the quarter cord, according to 


a = a m + aosin(u>t) (6) 

where a is the instantaneous angle of attack, a m is the mean angle of attack and a 0 is the pitch amplitude. The pitch frequency 
Co may be written in terms of a reduced frequency k r as Co = kjAU^/c, where c is the airfoil chord and Uoo is the free-steam 
velocity. Case 3 is defined by the parameter set a m = 4.86°, a m = 2.44°, k r = 0.01547, = 0.6, Re £ = 4.8xl0 6 . 

Within FUN3D, the angular displacement is specified as 
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where k c fd is a reduced frequency and t is the nondimensional time given by 


(7) 


^ t&oo 

Lref 

Recall that L re f is the dimensional reference length that makes the Reynolds number appearing in the stress tensor and heat 
flux vector of Eq. 1 a Reynolds number per unit length in the computational grid. To insure that the motion in the computational 
model tracks the experimental model, the arguments of the sin functions in Eqs. 6 and 7 must be identical. This leads to the 
value of k c fd in Eq. 7 
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The mesh used for this study was obtained by converting a structured mesh into an unstructured hexahedral mesh. In the 
mesh, the airfoil had a chord of unit length, so that L re f = c in Eq. 8. The original mesh contained 609 points around the 
airfoil and wake cut, and 257 points away from the airfoil. To evaluate the effect of grid refinement, this fine mesh was also 
coarsened successively in each direction to create meshes of size 305 x 129 and 153 x 65. For these simulations, FUN3D was 
run in 2D mode, in which the equations are solved on a single y=constant plane. This mode allows faster execution time by 
avoiding many unnecessary evaluations when there is no variation in the y direction. To assess the effect of time step on the 
solution, several different time steps were used, ranging from 100 time steps per pitch cycle to 800 time steps per pitch cycle. 
All computations of the 0012 airfoil were performed using the BDF 2 opt scheme and utilized the temporal error controller with 
a cutoff of 0.1. 


Figure 1 shows the computed lift coefficient vs angle of attack from the computations on the finest mesh, together with the 
experimental values from Ref. 15. The agreement is generally good, although there are relatively few experimental data points 
to determine whether the crossover seen in the computed results is seen in the experiment. While some variation is seen in the 
computation as the time step is refined from 100 to 200 steps per cycle, very little change is seen with further refinement to 400 
and 800 steps per cycle, indicating the computation is essentially temporally converged for this mesh. 

Figure 2 shows the corresponding pitching moment results. Pitching moment is always more difficult to accurately capture 
than lift coefficient, so it is not surprising to see more differences with the data than with lift coefficient. The largest discrepancy 
between computation and experiment occurs as the airfoil is pitching down (upper section of the curve), with the largest 
differences at the lower limit near a = 2.43°. In fact, as the time step is refined the computations move further away from 
the data in the pitch down region. However, the changes become smaller with decreasing time step, again indicating the 
computation is essentially temporally converged for this mesh. Very little variation with time step is seen in the pitch up region 
(lower section), for any time step. 

In Fig. 3 the effect of grid refinement on the lift variation is shown. For each mesh density, the smallest time step (800 
steps per cycle) was used. More variation in the result is seen with mesh refinement than with time step refinement. Figure 
4 illustrates the effect of grid refinement on pitching moment variation. Not unexpectedly, pitching moment is quite different 
between the coarse and medium grids. The change in pitching moment as the grid is refined from 305x129 to 609x257 is 
minimal, suggesting that an essentially mesh-size invariant solution is obtained for the finest mesh when using 800 steps per 
cycle. The fact that the computed moment curve on the finest mesh using the finest time step does not exactly replicate the 
experimental moment curve is perhaps not surprising, since the computation is strictly two dimensional, while no experiment 
can completely eliminate three-dimensional effects. 

Pressure coefficients at two points during the pitch cycle were extracted for comparison with the experimental data. The 
computed pressure data were extracted from the simulation corresponding to 800 time steps per cycle on the 609x257 mesh. 
Figure 5 shows pressure coefficients over the airfoil upper and lower surfaces as the airfoil is pitching up, at a = 5.92°. 
The computed and measured pressure coefficients are in excellent agreement. Figure 6 shows the comparison as the airfoil 
is pitching down, at a = 2.43°. At this point in the cycle, the agreement is generally very good, except for a slight under 
prediction of the peak suction as compared to the experiment. 


Pitching Blended Wing Body 

Next, a three-dimensional geometry is considered. The Blended Wing Body (BWB) is an advanced concept configuration 
that should offer an increase in F/D through reduced interference drag and a lower wetted area to volume ratio. 16 Several 
variants of the BWB have been wind-tunnel tested; the configuration considered here includes winglets and three nacelle/pylon 
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installations. Experimental measurements 17 were obtained with a relatively large model (approximately 7 foot wingspan) in 
the NASA Langley 14x22 subsonic tunnel. The model, with unpowered, flow-through nacelles, was designed for low speed 
measurements of static and dynamic stability derivatives. Data were taken for a number of flap/slat configurations; only the 
“clean” (no flaps/slats) geometry is considered here. As with any wind tunnel testing, mounting interference can be important. 
For the dynamic derivatives, the model was attached to a rather large mounting device that allowed the model to be pitched up 
and down in a sinusoidal fashion. The wind-tunnel model mounted on the pitch rig is shown in Figure 7. The top-mounted 
configuration is shown in the figure, however, the the experimental pitch data presented in this paper was taken with the model 
mounted from the bottom. Neither the mounting device nor the tunnel are modeled in the computation. Companion static force 
and moment computations, 18 including comparison with static data from another experimental facility, suggest the mounting 
system in the 14x22 tunnel imposed significant interference effects on the measured data. 

As with the NACA 0012, the pitching motion used in the BWB tests are described by Eq. 6. For the BWB, a range of a m , 
ao, and u> was tested in the wind tunnel. For the case presented here, a m = 0°, ao = 5°, and k r = Coc/2U 00 = 0.07, where 
c is the mean aerodynamic chord. Due to structural limitations, the tunnel dynamic pressure during the forced pitching had 
to be kept very low, around 2 psf, so that the resulting tunnel Mach number was approximately 0.036, with a corresponding 
.Reg = 0.7x10 6 . For the computational simulation, the Reynolds number was matched to the experimental value, but the Mach 
number was artificially boosted to 0.1 in order to avoid slow convergence of the compressible equations at very low Mach 
numbers. The sensitivity of the results to the change in Mach number from 0.036 to 0.1 was assumed to be negligible. It is 
important to note that this artificially larger Mach number must also be used when determining the reduced frequency k c fd 
from Eq. 8. Otherwise the ratio of characteristic time of the mean flow (e.g. the time it takes a fluid particle to traverse the 
body) and the characteristic time of the forced oscillation will be different between computation and experiment. 

The grid for this case was generated within the GridEx 19 framework. The mesh was developed in full-scale meters, so that 
so that L re f = 1 meter and c = 27.13228 meters. The initial mesh contains approximately 2.3 million nodes combined into 
tetrahedral cells. Source strengths were chosen to give a normal spacing near the wall leading to y + ~ 1 at the target Reynolds 
number. The results presented in this paper are for zero yaw angle, so only one half of the configuration is used, with symmetry 
boundary conditions on the center plane. Figure 8 shows the surface geometry and symmetry plane mesh. 

Figure 9 shows the pitching moment variation with angle of attack from both the computation and the experiment. Because 
the data is proprietary, labels are omitted from the vertical axes. Arrows labeled “+q" and “-q” indicate pitch up and pitch down 
segments of the curve, respectively. There is a clear offset between the data and the computation. This is very likely due to 
the significant post interference mentioned previously. The severity of this interference is indicated in the figure by symbols 
denoting the experimentally measured static pitching moment for a top-mounted model and a bottom-mounted model at the 
mean angle of attack. The dynamic data were taken in the wind tunnel with model mounted on the bottom. In the computational 
study, three time steps were used corresponding to 400, 800 and 1600 steps per pitch cycle. These results were all obtained 
using the second-order temporal scheme with an error controller cutoff of 0.1. Some variation of pitching moment with time 
step is seen between 400 and 800 steps per cycle, particularly in the range of a = +2° to a = +5°. However, very little 
change is observed in the pitching moment curve for time steps corresponding to 800 and 1600 steps per cycle, indicating that 
the pitching moment is essentially temporally converged for this mesh with 1600 steps per cycle. Apart from the static offset in 
the computed and measured pitching moment curves, the overall shape of the curve is in reasonably good agreement between 
the experiment and computation. 

Figure 10 shows the variation of normal force coefficient with angle of attack. Again there is a static offset from the 
experimental data, but overall the shape of the computational curve compares very favorably with the shape of the experimental 
curve. There is very little change in the normal force coefficient with time step; the largest change with time step refinement 
occurs in the range a = +2° toa = +5°, as was observed for the pitching moment coefficient. 

To assess the effect of mesh refinement on the solution, a second mesh was generated, containing 4.9 million nodes. The 
mesh was generated as an essentially uniform refinement of the initial grid. The variation of pitching moment coefficient with 
mesh size is shown in Fig. 11. Although the static offset from the experiment is slightly larger for the finer mesh, the overall 
shape of the pitching moment curve is somewhat better captured with the finer mesh. The change in solution between the 
two meshes is large enough that one cannot infer that the solution is mesh converged. Figure 12 shows the effect of mesh 
refinement on normal force coefficient, which tends to better align the slopes of the computational and experimental curves, 
without perceptible alteration to the shape of the computational curve. 

Next, computations for a m = 8° and a m = 16° were carried out on the 4.9 million node mesh, using the time step 
corresponding to 1600 steps per cycle. It should be noted that the range of Cm is the same in all the plots, to better indicate the 
relative variation; likewise the range of Cty is the same in all the plots. Figure 13 and 14 show the pitching moment and normal 
force coefficient variation, respectively, for a mean angle of attack of 8°. Again, the static values are also indicated in the plot 
to show the inherent offset between the two sets of experimental data and the computation. However, apart from the offset, the 
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computed pitching moment variation with angle of attack is in excellent agreement with the experimental data. The same is 
true for the normal force coefficient variation. The corresponding results for a mean angle of attack of a m = 16° are shown 
in Figs. 15 and 16. At this higher mean angle, the computed pitching moment variation is not particularly well captured by 
the computation, though the normal force coefficient variation, apart from the static offset, is still in reasonable agreement with 
the experimental data. Why the pitching moment curve does not show better agreement with the data is not certain, although 
flow visualization of the computed solutions show significantly more separation, and of a different nature at a m = 16° than at 
a m = 8°. It may well be that more grid points are required to adequately resolve the more complex flow field at a m = 16° 
before better agreement with the experimental data can be obtained. 

In stability and control applications, the so-called stability derivatives are widely used. These derivatives arise out of the 
assumption of small deviations from a reference condition of steady flight. Because of the small perturbation assumption, 
stability derivatives are typically valid over a limited region around the reference condition; in highly non-linear flight regimes 
they may not be valid at all. Because these derivatives are so widely used, and because in many cases they are difficult to obtain 
experimentally, it is desirable to use computational methods to estimate them. For longitudinal dynamic stability assessment, a 
parameter of great interest is the “lumped pitch damping derivative” 

j, dC M 

Cm «-^T 

where q = qc/{2U 00 ). The overbar indicates that this derivative contains both pitch rate (q) and angle of attack rate (d) 
terms. In a wind tunnel or flight test (or computation) in which the configuration undergoes pitching about a spanwise axis 
(with associated angle of attack change), it is not possible to separate the q effect from the d effect. In fact, the separation of 
the two effects is another consequence of the linearized, small perturbation assumption. For a stable aircraft. Cm? should be 
negative, that is, a positive pitch rate (pitch up) should generate a nose-down moment. The most straightforward method for 
determining Cm q is to use a finite -difference between Cm values as measured (or computed) while the model is pitching up and 
pitching down, at the reference (mean) angle of attack. Table 2 shows the relative errors between C’Mq found by applying finite 
differences to both the experimental and computational Cm data, at the three mean angles of attack considered above. Owing 
to the proprietary nature of the data, percent errors, relative to the experimental data, are given. At each value of a m , CMq was 
found to be negative, indicating a dynamically stable vehicle in pitch over this range of angle of attack. It is interesting to note 
that although the overall agreement between the computed and measured moment curves is better at a m = 8° than at a m = 0° 
(and significantly better than at a m = 16°), the computed value of C Mq at a m = 8° is much further away from the measured 
value than at a m = 0°, and only in somewhat better agreement than at a m = 16°. The excellent agreement between computed 
and measured C Mq at a m = 0° is almost certainly fortuitous, and illustrates the pitfall of using a single derived value as a 
measure of the quality of the computational result. 



% Error 

0° 

-4 

8° 

+37 

16° 

-45 


Table 2 Errors in computed values of Cm, 


Synthetic Jet 

The final application presented in this paper is a piston-driven synthetic jet. A schematic of the actuator used in the 
experiment is shown in Figure 17. Owing to the vertical movement of the piston, flow passed in and out of the cavity through 
an orifice. It should be noted that while the orifice was of circular cross section, the piston consisted of a rigid square plate 
attached somewhat asymmetrically to the flexible membrane. Thus the flow into and out of the orifice was not completely 
symmetric about the centerline. The orifice was located on a splitter plate mounted in a wind tunnel. The Mach number in the 
tunnel was 0.1, and the boundary layer height near the orifice of approximately 21mm. It should be noted that the cavity was 
quite shallow relative to its width, and the piston movement caused a considerable change in volume as it cycled up and down. 

The geometry of the cavity is simplified by modeling the piston/membrane as a flat surface, with the vertical position of the 
surface taken such that the “at rest” volume of the cavity is unchanged. To reduce the problem size, the domain is divided at 
the tunnel centerline and only half of the domain is modeled in the computation. The computational domain extends from 8 jet 
diameters upstream and to the side of jet center, and 16 jet diameters downstream of the jet center. 
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The grid for this case was also generated using GridEx, and the wall spacing set so that the resulting y + ~ 1 for the target 
conditions. Grid spacings inside the cavity were based on preliminary simulations with steady blowing applied at the piston 
face. Griding sources were placed around the lip of the jet and along the jet centerline in an effort to improve the clustering in 
these regions. The resulting mesh contained approximately 255,000 nodes. Fig. 18 shows the grid in symmetry plane for this 
configuration. To assess the effect of mesh size on the computed solutions, a second mesh of approximately 46,000 nodes was 
constructed as an essentially uniform coarsening of the finer mesh. 

In order to match the conditions in the tunnel, FUN3D was modified to allow a velocity profile to be input. A turbulent 
velocity profile obtained with FUN3D and the SA turbulence model on a flat plate grid with conditions that gave a thickness of 
approximately 21mm was used as an inflow condition. 

Rather than modeling the piston motion, at the mean piston surface a time-dependent velocity was imposed, given in 
dimensional form by 


u = 0 0 = 0 W = [(pw) max / p\cOS (2n ft) 

where / is the piston frequency, 150 Hz. Density and pressure on this boundary were left to evolve as part of the solution 
process. The value of (pt0) max is chosen as part of a calibration process to approximately match the peak vertical velocity at the 
outflow of the orifice between computation and experiment. The value of this calibration constant is (pw) max = O.OOTpooCtoo. 
All results were obtained with a time step corresponding to 720 steps per actuator cycle. 

The experimental data for this configuration at several locations, marked on Fig. 19, are available from the experimental 
study of Schaeffer and Jenkins. 20 Figures 20 and 21 show the comparison of the computed phase-averaged normal and stream- 
wise velocities with the experimental data at .7=50.8 and 2=0.4 mm location, which corresponds to the center of the jet exit. The 
overall agreement with the experimental data is fair and present results are comparable to numerical results of other researchers 
who employed similar boundary conditions 4 to simulate this case. Very little change with mesh refinement is observed in the 
normal velocity with mesh refinement although the streamwise component shows some improvement with the finer mesh for 
phase angles between approximately 80° and 220°. 

The phase-averaged velocities at a distance of 2 jet diameters downstream of the jet centerline ( 7=63.5 mm) are compared 
with the experimental data in Figs. 22 and 23. The computational results on the finer mesh compare favorably with the 
experimental data at this location; the features on the coarser mesh are greatly dissipated. Finally, the time-averaged velocity 
profiles at 7=63.5 mm are shown in Figs. 24 and 25. The computational results replicate the experimentally observed features 
and are again comparable to the CFD results of other researchers. 4 The time-averaged normal velocity shows the greatest 
sensitivity to mesh refinement. The peak velocity is under predicted on the coarser mesh and over predicted on the finer mesh, 
although the profile on the finer mesh is in better agreement overall with the experimental data. 

Summary 

A Navier-Stokes solver for unstructured grids has been extended to allow time-accurate computations on grids that move 
as a rigid body. The time-accurate capability of the solver has been extended to third-order accurate BDF schemes, including 
a temporal-error controller that eliminates much of the guesswork involved in determining the proper number of nonlinear 
subiterations required within a given time step. In addition, a transpiration boundary condition for simplified modeling of 
synthetic jet flows has been developed. 

The modified code has been applied to several test cases of interest to the active flow control and stability and control 
communities. The influence of time step and mesh size on the solutions was assessed. These initial results show reasonable 
agreement with available experimental data in most cases. Both the BWB and synthetic jet cases presented here involved 
simplifications of the physical model - neglect of mounting hardware or tunnel walls in the BWB case and use of a transpiration 
boundary condition rather than a moving piston in the synthetic jet case. Recently the mesh movement within the FUN3D solver 
has been extended to include deforming meshes, and the computations for these cases will be repeated taking the mounting 
system into account for the BWB and the moving piston into account for the synthetic jet. 
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Fig. 1 Effect of time-step refinement on iift coefficient vs. Fig. 2 Effect of time-step refinement on pitching moment 

angle of attack, 0012 airfoil, AGARD Case 3. coefficient vs. angle of attack, 0012 airfoil, AGARD Case 3. 
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Fig. 3 Effect of grid refinement on lift coefficient vs. angle Fig. 4 Effect of grid refinement on pitching moment coeffi- 

of attack, 0012 airfoil, AGARD Case 3. cient vs. angle of attack, 0012 airfoil, AGARD Case 3. 
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Fig. 5 Pressure coefficient distribution during pitch-up, 
a = 5.92°, 0012 Airfoil, AGARD Case 3. 



Fig. 6 Pressure coefficient distribution during pitch-down, 
a = 2.43°, 0012 Airfoil, AGARD Case 3. 



Fig. 7 3% scale model Blended Wing Body on the dynamic Fig. 8 Computational model of Blended Wing Body; 2.3 

test rig in the 14 x 22 Wind Tunnel; top mounting shown. million node mesh. 
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a, deg. a, deg. 


Fig. 9 Effect of time-step refinement on BWB pitching mo- Fig. 10 Effect of time-step refinement on BWB normal force 

ment coefficient vs. angle of attack, a m = 0°; “+q” and “-q” coefficient vs. angle of attack, a m = 0°. 

indicate the pitch-up and pitch-down legs of the oscillation. 




a, deg. a, deg. 


Fig. 11 Effect of grid refinement on BWB pitching moment 
coefficient vs. angle of attack, a m = 0°. 


Fig. 12 Effect of grid refinement on BWB normal force co- 
efficient vs. angle of attack, a m = 0°. 
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a, deg. a, deg. 


Fig. 13 BWB pitching moment coefficient vs. angle of attack, 

a m = 8 °. 



a, deg. 


Fig. 14 BWB normal force coefficient vs. angle of attack, 
= 8 °. 



a, deg. 


Fig. 15 BWB pitching moment coefficient vs. angle of attack, Fig. 16 BWB normal force coefficient vs. angle of attack, 

a m = 16°. a m = 16°. 
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Fig. 17 Sketch of fluidic actuator, showing cavity and piston (from Ref. 21). 
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Fig. 18 Symmetry plane grid for the actuator. Fig. 19 Sketch showing locations of data comparison with 

experiment (from Ref. 21). 
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Fig. 20 Normal velocity component as a function of actuator 
phase angle, x=50.8 mm. 


Fig. 21 Streamwise velocity component as a function of ac- 
tuator phase angle, x=50.8 mm. 




Fig. 22 Normal velocity component as a function of actuator Fig. 23 Streamwise velocity component as a function of ac- 

phase angle, x=63.5 mm. tuator phase angle, x=63.5 mm. 
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mm 




Fig. 24 Time averaged normal velocity component, x=63.5 Fig. 25 Time averaged streamwise velocity component, 

mm. x=63.5 mm. 
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